#include "Data.hh"

using namespace std;

//ClassImp(Element);
//ClassImp(Material);
//ClassImp(Isotopes);

Data::Data()
{
}

Data::Data(int argc, char* argv[])
{
  settings = NULL;
  HCName = NULL;
  IsotopeTable = NULL;
  TargetMaterial = NULL;
  ElasticTarget = NULL;
  LevelEnergy = NULL;
  LevelWidth = NULL;
  FeedingProbability = NULL;
  Factor = NULL;
  Exponent = NULL;
  NumberOfDecays = NULL;
  LevelSpin = NULL;
  LevelParity = NULL;
  LevelID = NULL;
  DecayLevel = NULL;
  DecayProbability = NULL;
  DecayType = NULL;
  DecayDelta = NULL;

  TargetRatio = 0.1;

  int i;

  //setting directories
  SimulationDirectory = "/usr/local/TREX_simulation";
  MassFile = SimulationDirectory;
  MassFile.append("/cmass.dat");

  ScratchDirectory = "/Users/JacobSJ/Physics/is430/2010/Simulation";

  //set default values
  ISession = 0;
  Vis      = 0;
  NoE      = 1000;
  Filename = "test.root";
  BeamWidth = 5.*mm;
  AngleWidth = 0.*degree;
  TestFlag = 0;
  Tracking = 0;
  CDhit = 0;
  BeamEnergy = 90.*MeV;
  Settingsfile = SimulationDirectory+"/Settings/DefaultSettings.dat";
  TrackAll = 0;
  WriteAll = false;
  BarrelDetectors = 4;
  Quadrants = 4;
  BarrelType = 1;
  CDPad = 0;
  NumberOfLevel = 0;
  LevelFile = SimulationDirectory+"/LevelFiles/DefaultLevels.dat";
  BetaFile = SimulationDirectory+"/BetaLevels/DefaultLevels.dat";
  OpticalPotentialFile = "";
  OpticalPotential = false;
  ThicknessStepSize = 0;
  ProjectileName = "30Mg";
  TargetName = "2H";
  TransferFlagN = 1;
  TransferFlagP = 0;
  NumberOfThicknessSteps = 1024;
  CDResolution = 50.;//*keV;
  BarrelResolution = 50.;//*keV;
  TransferCrossSection = 1;
  Cut = 0;
  ThetaCmMin = 0.*degree;
  ThetaCmMax = 180.*degree;
  DetectionEnergy = 0.*MeV;
  EuropiumSource = false;
  CobaltSource = false;
  ChamberType = 4;
  Gammas = true;
  Efficiency = false;
  Coulex = false;
  IsoGamma = false;
  DownScale = 0;
  OnePartOnly = 0;
  for(i=0;i<8;i++)
    {
      GammaCounter[i] = 0;
    }
  VerboseLevel = 0;
  IsotropeParticleDistribution = false;
  DeltaE = false;
  RotateMiniball = false;
  ExperimentType = 0;
  Screening = false;
  ScreeningMaterialName.erase();
  AlphaSource = false;
  QuadAlphaSource = false;
  BetaEnergy = 0.0;
  AddBeta = 0.0;
  RatioBeta = 0.0;
  BackwardFoil = 0;
  noCarbon = 0;
  Lifetime = -1;
  G4String tmp;

  for(i=0;i<argc;i++)
    {
      cout<<argv[i]<<" ";
      CommandLine += argv[i];
      CommandLine += " ";
    }
  cout<<endl;

  vector<string> reaction_tmp;
  vector<double> addbeta_tmp;
 
  //instance of command line interface used to read the command line flags 
  CommandLineInterface* interface = new CommandLineInterface();

  interface->Add("-a",  "write all events, detected or not", &WriteAll);
  interface->Add("-ab", "add beta source, i.e. transfer and beta source <Q value [MeV]> <ratio beta/transfer>", &addbeta_tmp);  
  interface->Add("-aq", "simulate quadruple alpha source with 148Gd, 239Pu, 241Am and 244Cm", &QuadAlphaSource);
  interface->Add("-as", "simulate triple alpha source with 239Pu, 241Am and 244Cm", &AlphaSource);
  interface->Add("-be", "<beam energy>                          (default is 90 MeV)", &BeamEnergy);
  interface->Add("-bf", "<backward barrel foil>                 (0 = no foil, 1 = cylinder, 2 = perpendicular to beam, 3 = perpendicular to beam + backward CD)", &BackwardFoil);
  interface->Add("-bl", "level file for betas", &BetaFile);
  interface->Add("-br", "<barrel resolution>                    (fwhm in keV)", &BarrelResolution);
  interface->Add("-bs", "simulate beta source, <Q value [MeV]>", &BetaEnergy);
  interface->Add("-bt", "barrel type, 0 = 2007 setup, no BE detectors\n            1 = 2008 setup with BE detectors", &BarrelType);
  interface->Add("-bw", "<beam width>                           (set FWHM of the beam in mm)", &BeamWidth);
  interface->Add("-Aw", "<angle width>                           (set FWHM of the beam in degree)", &AngleWidth);
  interface->Add("-c",  "<cut type>                             (select different types of cuts:\n                                         1 = chamber is cut open,                      2 = no Miniball except behind chamber,\n                                         4 = no left barrel detectors,        8 = no left segment of CDs,\n                                        16 = no left segment of PCB of CDs,           32 = no Barrel PCB,\n                                       and combinations of above)", &Cut);
  interface->Add("-cd", "only track if hit in CD", &CDhit);
  interface->Add("-cl", "Coulex", &Coulex);
  interface->Add("-co", "simulate 60Co source instead of particles", &CobaltSource);
  interface->Add("-cp", "0 no CDE Detectors, 1 CDE detectors", &CDPad);
  interface->Add("-cr", "<CD resolution>                        (fwhm in keV)", &CDResolution);
  interface->Add("-ct", "<chamber type>                         (set ChamberType:\n                                       0 = no Chamber,                               1 = small cylindrical Chamber,\n                                       2 = middle cylindrical Chamber,               3 = long cylindrical Chamber,\n                                       4 = spherical Chamber,                        5 = new cylindrical Chamber,\n                                      10 = old Chamber)", &ChamberType);
  interface->Add("-de", "<detection energy>                     (minimum energy for event to be written away in MeV)", &DetectionEnergy);
  interface->Add("-ds", "down scaling with 2^(factor)", &DownScale);
  interface->Add("-dt", "Test of DeltaE, with isotropic lab distribution", &DeltaE);
  interface->Add("-ef", "simulate gammas for efficiency", &Efficiency);
  interface->Add("-eu", "simulate 152 Eu source instead of particles", &EuropiumSource);
  interface->Add("-ex", "<experiment_type>                      (set Eperiment Type:\n                                       0 = one CD forward,                          1 = Mahmoud's CD\n                                       2 = only barrel,                              3 = barrel + backward CD\n                                       4 = barrel + 2 CDs", &ExperimentType);
  interface->Add("-f",  "<filename>", &Filename);
  interface->Add("-g",  "generate gammas", &Gammas);
  interface->Add("-i",  "interactiv", &ISession);
  interface->Add("-ig", "isotrope gamma distribution", &IsoGamma);
  interface->Add("-ip", "isotrope particle distribution", &IsotropeParticleDistribution);
  interface->Add("-lf", "level_file", &LevelFile);
  interface->Add("-lt", "lifetime in ns", &Lifetime);
  interface->Add("-n",  "<# of events>", &NoE);
  interface->Add("-nc",  "scattering on carbon 0 yes, 1 no", &noCarbon);
  interface->Add("-op", "<optical potential file>               (file with optical potential for all target elements)", &OpticalPotentialFile);
  interface->Add("-pt",  "simulate only 1 ejectile, 2 recoil, 0 all", &OnePartOnly);  
  interface->Add("-r",  "<projectile> <target> <# of neutrons>  (reaction of projectile with target,\n                                        projectile gains # of neutrons in reaction)", &reaction_tmp);
  interface->Add("-rm", "rotate miniball angles", &RotateMiniball);
  interface->Add("-sc", "<Material>                             (screening of detectors with 'Material' foil)", &ScreeningMaterialName);
  interface->Add("-sf", "<settings file>                        (specify file with the settings)", &Settingsfile);
  interface->Add("-t",  "run test part of simulation", &TestFlag);
  interface->Add("-tc", "<transfer cross section>               (negative -> calculate, otherwise 0->1)", &TransferCrossSection);
  interface->Add("-tm", "<theta_cm min.>                        (min. value of theta_cm in degree)", &ThetaCmMin);
  interface->Add("-tM", "<theta_cm max.>                        (max. value of theta_cm in degree)", &ThetaCmMax);
  interface->Add("-ts", "<# of ThicknessSteps>                  (set number of steps used in EVsX and RVsX)", &NumberOfThicknessSteps);
  interface->Add("-v",  "turn on visualization", &Vis);
  interface->Add("-vl", "<verbose level>", &VerboseLevel);

  if(!interface->CheckFlags(argc, argv)){
    exit(1);
  }
  cout << "-----------------------------------------------------"<<endl;
  cout << "lifetime " << Lifetime << " ns "  << endl;
  cout << "-----------------------------------------------------"<<endl;
  
  //if(Lifetime>-0.5){
  //  Efficiency=true;
  //}

  if(argc == 1)
    {
      exit(1);
    }
  if(!Filename.contains(".root"))
    {
      Filename += ".root";
    }

  if(!OpticalPotentialFile.empty())
    {
      OpticalPotential = true;
    }

  if(!ScreeningMaterialName.empty())
    {
      Screening = true;
      ScreeningMaterialName.toLower();
    }

  if(Vis)
    {
      ISession = true;
    }
  if(noCarbon>1){
    exit(1);
  }


  if(reaction_tmp.size() == 3)
    {
		ProjectileName = reaction_tmp[0];
      TargetName =     reaction_tmp[1];
      
		if(atoi(reaction_tmp[2].c_str())>=10)
		{
			TransferFlagN = -1*atoi(reaction_tmp[2].c_str())/10;
			TransferFlagP = -1*atoi(reaction_tmp[2].c_str())%10;
			cout << "HER::: "<< TransferFlagN << "	" << TransferFlagP << endl;
		}
		else
		{
			TransferFlagN =   atoi(reaction_tmp[2].c_str());
    	}
	 }
  else if(reaction_tmp.size() != 0)
    {
      cerr<<"flag -r provided but not three arguments following: ";
      for(i = 0; i < reaction_tmp.size(); i++)
	{
	  cerr<<reaction_tmp[i]<<" ";
	}
      exit(1);
    }
  if(addbeta_tmp.size() == 2)
    {
      AddBeta = addbeta_tmp[0];
      RatioBeta = addbeta_tmp[1];
    }
  else if(addbeta_tmp.size() != 0)
    {
      cerr<<"flag -ab provided but not two arguments following: ";
      for(i = 0; i < addbeta_tmp.size(); i++)
	{
	  cerr<<addbeta_tmp[i]<<" ";
	}
      exit(1);
    }
  

  if(VerboseLevel > 2)
    {
      cout<<"ReadArguments finished without error"<<endl;
    }
  if(DownScale>0&&Coulex){
    cout << "DownScale = " << DownScale << " supressing no gamma events with " << pow(2.,(int)DownScale) << endl;
    DownScale = pow(2.,(int)DownScale);
  }
  else
    DownScale=0;

  if(CobaltSource && EuropiumSource)
    {
      cout<<"Simulate cobalt source AND europium source? I'm confused!"<<endl;
      exit(3);
    }
  if(BackwardFoil && !Screening)
    {
      cout<<"Foil geometry given, but screening deactivated??? Set -bf to 0 or choose screening material -sc"<<endl;
      exit(3);
    }
  if(BarrelType == 1)
    BarrelDetectors = 4;
  else if(BarrelType == 0)
    BarrelDetectors = 3;
  else{
    cerr << "bt must be 0 or 1" << endl;
    exit(3);
  }

  if(ExperimentType == 0){
    FirstBarrelHC = 4;
    BarrelDetectors = 0;
  }
  if(ExperimentType == 1){
    FirstBarrelHC = 4;
    BarrelDetectors = 0;
  }
  else if(ExperimentType == 2){ //no CDs
    FirstBarrelHC = 0;
  }
  else if(ExperimentType == 3){ //one CD + barrel 
    FirstBarrelHC = 4;
  }
  else if(ExperimentType == 4){ //FCD,BCD, +4 F,B,FE, BE
    FirstBarrelHC = 8;
    if(CDPad)
      FirstBarrelHC+=4;
  }
  if(CDPad)
    FirstBarrelHC+=4;

  if(ChamberType == 10)
    {
      tmp = SimulationDirectory+"/Settings/DefaultSettings.dat";
      if(Settingsfile == tmp.data())
	{
	  Settingsfile = SimulationDirectory+"/Settings/OldChamberSettings.dat";
	}
      else
	{
	  cout<<"WARNING: You chose to simulate the old setup, but you've also demanded to load settings from "<<Settingsfile<<endl
	      <<"I'm ignoring those settings and will read them from "<<SimulationDirectory<<"/OldChamberSettings.dat instead"<<endl;
	  //Settingsfile = SimulationDirectory+"/Settings/OldChamberSettings.dat";
	}
      if(BarrelDetectors!=0)
	{
	  cout<<"WARNING: You chose to simulate the old setup, but you've also demanded to simulate "<<BarrelDetectors<<" barrel detectors"<<endl
	      <<"I'm ignoring that and will simulate no barrel detectors"<<endl;
	  BarrelDetectors = 0;	
	}
      FirstBarrelHC = 4;
    }


  NumberOfHC = Quadrants*BarrelDetectors+FirstBarrelHC;

  cout << "FirstBarrelHC " << FirstBarrelHC << " NumberOfHC " << NumberOfHC << endl;

  settings = new Settings(Settingsfile,VerboseLevel,Cut,RotateMiniball);

  HCName = new G4String[NumberOfHC];

  if(ReadLevels() != 0){
    cerr<<"Error in Data::ReadLevels()!"<<endl;
    exit(1);
  }
  if(ReadBetaLevels() != 0){
    cerr<<"Error in Data::ReadBetaLevels()!"<<endl;
    exit(1);
  }

  if(VerboseLevel > 0)
    {
      cout<<"read level file "<<LevelFile<<" with "<<NumberOfLevel<<" levels, AngularDistributionFile = "<<AngularDistributionFile<<":"<<endl;
      for(i=0;i<NumberOfLevel;i++)
	{
	  cout<<"\t"<<LevelEnergy[i]/keV<<" "<<FeedingProbability[i]<<" "<<Factor[i]<<" "<<Exponent[i]<<endl;
	}

      cout<<"using mass file "<<MassFile<<endl;
    }

  IsotopeTable = new Isotopes(MassFile.data());

  Projectile = IsotopeTable->Search((char*)ProjectileName.data());
  Target = IsotopeTable->Search((char*)TargetName.data());
  if(Coulex){
    Ejectile = IsotopeTable->Search((char*)ProjectileName.data());
    Recoil = IsotopeTable->Search((char*)TargetName.data());
    TransferCrossSection = 1;
  }
  else{
    Ejectile = IsotopeTable->Reaction(Projectile,TransferFlagP,TransferFlagN);
    Recoil = IsotopeTable->Reaction(Target,-TransferFlagP,-TransferFlagN);
  }
  if(Projectile == NULL || Ejectile == NULL || Target == NULL || Recoil == NULL)
    {
      cout<<"failed to find Isotopes in IsotopeTable ("<<MassFile<<"):"<<Projectile<<" "<<Ejectile<<" "<<Target<<" "<<Recoil<<endl;

      exit(1);
    }

  //PE and MY are implemented as materials, everything else should be the name of the element
  if(settings->TargetMaterialName.contains("PE") || settings->TargetMaterialName.contains("MY"))
    {
      TargetMaterial = new Material((char*)settings->TargetMaterialName.data());
    }
  else
    {
      //if target material name is the same as the name of the scattering target build set the material to only this element
      if(settings->TargetMaterialName == TargetName)
	{
	  TargetMaterial = new Material((char*)settings->TargetMaterialName.data(),false);
	}
      else
	{
	  if(Coulex){
	    cerr << " Error, Coulex only possible with one isotopes in target " << endl;
	    exit(3);
	  }
	  char* ElementNames[] = {(char*)settings->TargetMaterialName.data(),(char*)TargetName.data()};
	  double ElementRatios[] = {1-TargetRatio,TargetRatio};
	  TargetMaterial = new Material(2,ElementNames,ElementRatios,false);
	}
    }

  ElasticTarget = new Element*[TargetMaterial->NumberOfElements()];
  for(i=0;i<TargetMaterial->NumberOfElements();i++)
    {
      ElasticTarget[i] = TargetMaterial->GetElement(i);
    }

//  if(TransferCrossSection < 0)
//    {
//      cout<<"Warning: calculating TransferCrossSection is not implemented yet (setting to 1)!"<<endl;
//      TransferCrossSection = 1;
//    }
  
  if(NumberOfThicknessSteps < 16)
    {
      NumberOfThicknessSteps = 16;
      if(VerboseLevel > 3)
	cout<<"changed NumberOfThicknessSteps to "<<NumberOfThicknessSteps<<endl;
    }
  ThicknessStepSize = settings->TargetThickness/NumberOfThicknessSteps;
  
  EjectileParticleName = Ejectile->Name() + G4UIcommand::ConvertToString(Ejectile->A());
  RecoilParticleName = Recoil->Name() + G4UIcommand::ConvertToString(Recoil->A());


}

Data::~Data()
{
  if(settings != NULL)
    delete settings;
  if(HCName != NULL)
    delete[] HCName;
  if(IsotopeTable != NULL)
    delete IsotopeTable;
  if(TargetMaterial != NULL)
    delete TargetMaterial;
  if(ElasticTarget != NULL)
    delete[] ElasticTarget;
  if(LevelEnergy != NULL)
    delete[] LevelEnergy;
  if(LevelWidth != NULL)
    delete[] LevelWidth;
  if(FeedingProbability != NULL)
    delete[] FeedingProbability;
  if(BetaLevelEnergy != NULL)
    delete[] BetaLevelEnergy;
  if(BetaPopulation != NULL)
    delete[] BetaPopulation;
  if(Factor != NULL)
    delete[] Factor;
  if(Exponent != NULL)
    delete[] Exponent;
  if(NumberOfDecays != NULL)
    delete[] NumberOfDecays;
  if(LevelSpin != NULL)
    delete[] LevelSpin;
  if(LevelParity != NULL)
    delete[] LevelParity;
  if(LevelID != NULL)
    delete[] LevelID;
  for(int i=0;i<NumberOfLevel;i++)
    {
      if(DecayLevel != NULL)
	delete[] DecayLevel[i];
      if(DecayProbability != NULL)
	delete[] DecayProbability[i];
      if(DecayType != NULL)
	delete[] DecayType[i];
      if(DecayDelta != NULL)
	delete[] DecayDelta[i];
    }
  if(DecayLevel != NULL)
    delete[] DecayLevel;
  if(DecayProbability != NULL)
    delete[] DecayProbability;
  if(DecayType != NULL)
    delete[] DecayType;
  if(DecayDelta != NULL)
    delete[] DecayDelta;
}

int Data::ReadBetaLevels()
{
  int i;
  double pop_tot=0;

  cout<<"start reading "<<BetaFile<<" ... \n";
  ifstream file(BetaFile);
  if(file.bad())
    {
      cerr<<"Unable to open "<<BetaFile<<"!\nexiting ... \n";
      exit(2);
    }
  file.ignore(1000, '\n');
  file.ignore(1000, '\n');
  file>>BetaNrOfLevels;
  file.ignore(1000, '\n');
  if(VerboseLevel > 2){
      cout<<"reading in "<<BetaNrOfLevels<<" levels"<<endl;
    }
  file.ignore(1000, '\n');

  BetaLevelEnergy = new Double_t[BetaNrOfLevels];
  BetaPopulation = new Double_t[BetaNrOfLevels];
  
  for(i=0;i<BetaNrOfLevels;i++){
    file >> BetaLevelEnergy[i] >> BetaPopulation[i];
    if(VerboseLevel > 2){
      cout << "level " << i << ": E = " << BetaLevelEnergy[i] << ", population = " << BetaPopulation[i] << endl;
    }
    BetaLevelEnergy[i] *= keV;
    pop_tot += BetaPopulation[i];
    
    file.ignore(1000, '\n');
  }
  for(i=0;i<BetaNrOfLevels;i++){
    BetaPopulation[i] /= pop_tot;
    if(i>0)
      BetaPopulation[i] += BetaPopulation[i-1];
  }  

  file.close();

  return 0;

}

int Data::ReadLevels()
{
  int i,j;
  G4double tmp = 0;
  cout<<"start reading "<<LevelFile<<" ... \n";
  ifstream file(LevelFile);
  if(file.bad())
    {
      cerr<<"Unable to open "<<LevelFile<<"!\nexiting ... \n";
      exit(2);
    }
  file.ignore(1000, '\n');
  file.ignore(1000, '\n');
  file>>AngularDistributionFile;
  file.ignore(1000, '\n');
  file>>NumberOfLevel;
  file.ignore(1000, '\n');
  if(VerboseLevel > 2){
    cout<<"reading in "<<NumberOfLevel<<" levels"<<endl;
  }
  LevelEnergy = new Double_t[NumberOfLevel];
  LevelWidth = new Double_t[NumberOfLevel];
  FeedingProbability = new Double_t[NumberOfLevel];
  Factor = new G4double[NumberOfLevel];
  Exponent = new G4double[NumberOfLevel];
  if(Gammas)
    {
      NumberOfDecays = new G4int[NumberOfLevel];
      LevelSpin = new G4double[NumberOfLevel];
      LevelParity = new G4int[NumberOfLevel];
      LevelID = new G4int[NumberOfLevel];
      DecayLevel = new G4int*[NumberOfLevel];
      DecayProbability = new G4double*[NumberOfLevel];
      DecayType = new G4int*[NumberOfLevel];
      DecayDelta = new G4double*[NumberOfLevel];
    }
  for(i=0;i<NumberOfLevel;i++)
    {
      file>>LevelEnergy[i]>>LevelWidth[i]>>FeedingProbability[i]>>Factor[i]>>Exponent[i];
      if(VerboseLevel > 2){
	cout<<"level "<<i<<": E = "<<LevelEnergy[i]<<" ("<<LevelWidth[i]<<") "<<", feeding = "<<FeedingProbability[i]<<" ("<<Factor[i]<<", "<<Exponent[i]<<")"<<endl;
      }
      LevelEnergy[i] *= keV;
      LevelWidth[i] *= keV;
      tmp += FeedingProbability[i];
      if(Gammas)
	{
	  file>>LevelSpin[i]>>LevelParity[i]>>NumberOfDecays[i];
	  if(VerboseLevel > 2)
	    {
	      cout<<"adding "<<NumberOfDecays[i]<<" decays to level "<<i<<" (spin = "<<LevelSpin[i]<<", parity = "<<LevelParity[i]<<")"<<endl;
	    }
	  DecayLevel[i] = new G4int[NumberOfDecays[i]];
	  DecayProbability[i] = new G4double[NumberOfDecays[i]];
	  DecayType[i] = new G4int[NumberOfDecays[i]];
	  DecayDelta[i] = new G4double[NumberOfDecays[i]];
	  for(j=0;j<NumberOfDecays[i];j++)
	    {
	      file>>DecayLevel[i][j]>>DecayProbability[i][j]>>DecayType[i][j]>>DecayDelta[i][j];
	      if(VerboseLevel > 2)
		{
		  cout<<"read decay "<<j<<": decaying to level "<<DecayLevel[i][j]<<", with probability "<<DecayProbability[i][j]<<", type = "<<DecayType[i][j]<<", delta = "<<DecayDelta[i][j]<<endl;
		}
	    }
	}

      file.ignore(1000, '\n');
    }

  if(VerboseLevel > 1)
    {
      cout<<"read level file "<<LevelFile<<" with "<<NumberOfLevel<<" excitation energies and uncorrected probabilities (sum = "<<tmp<<"):"<<endl;
      cout<<"AngularDistributionFile = "<<AngularDistributionFile<<endl;
      for(i=0;i<NumberOfLevel;i++)
	{
	  cout<<"\t"<<LevelEnergy[i]/keV<<" "<<FeedingProbability[i]<<" "<<Factor[i]<<" "<<Exponent[i];
	  if(Gammas)
	    {
	      cout<<" "<<NumberOfDecays[i];
	      for(j=0;j<NumberOfDecays[i];j++)
		cout<<" "<<DecayLevel[i][j]<<" "<<DecayProbability[i][j]<<" "<<DecayType[i][j]<<" "<<DecayDelta[i][j];
	    }
	  cout<<endl;
	}
    }

  //normalize probability to sum = 1
  for(i=0;i<NumberOfLevel;i++)
    {
      FeedingProbability[i] /= tmp;
      if(i>0)
	FeedingProbability[i] += FeedingProbability[i-1];
    }

  file.close();

  return 0;
}

void Data::run(G4RunManager* runManager)
{
  if(ISession == 1)
    {
      // Define UI terminal for interactive mode
      G4UIsession* session = new G4UIterminal(new G4UItcsh);

      if(Vis == 1)
	{
	  G4UImanager* UI = G4UImanager::GetUIpointer();     
	  G4String tmp = "/control/execute "+SimulationDirectory+"/OGL_vis.mac";
	  UI->ApplyCommand(tmp.data()); 
	}

      session->SessionStart();
      delete session;
    }
  else if(NoE > 0 || TestFlag)
    {
      runManager->BeamOn(NoE);
    }
}
void Data::CreateBranches(TTree* tree)
{
  tree->Branch("BeamWidth",&BeamWidth,"BeamWidth/D");
  tree->Branch("BeamEnergy",&BeamEnergy,"BeamEnergy/D");
  tree->Branch("NumberOfBarrelDetectors",&BarrelDetectors,"NumberOfBarrelDetectors/I");
  tree->Branch("CDResolution",&CDResolution,"CDResolution/D");
  tree->Branch("BarrelResolution",&BarrelResolution,"BarrelResolution/D");
  tree->Branch("TransferCrossSection",&TransferCrossSection,"TransferCrossSection/D");
  tree->Branch("ThetaCmMin",&ThetaCmMin,"ThetaCmMin/D");
  tree->Branch("ThetaCmMax",&ThetaCmMax,"ThetaCmMax/D");
  tree->Branch("DetectionEnergy",&DetectionEnergy,"DetectionEnergy/D");
  tree->Branch("NumberOfLevel",&NumberOfLevel,"NumberOfLevel/I");
  char tmp[256];
  sprintf(tmp,"LevelEnergy[%d]/D",NumberOfLevel);
  tree->Branch("LevelEnergy",LevelEnergy,tmp);
  tree->Branch("NoE",&NoE,"NoE/I");
  tree->Branch("TransferFlag",&TransferFlagN,"TransferFlag/I");
  tree->Branch("ProjectileName",(char*) ProjectileName.data(),"ProjectileName/C");
  tree->Branch("TargetName",(char*)TargetName.data(),"TargetName/C");
  tree->Branch("TargetMaterialName",(char*)settings->TargetMaterialName.data(),"TargetMaterialName/C");
  settings->CreateBranches(tree);
}
int Data::DetectorType(int HCID){
  //returns:
  //0 forward barrel
  //1 backward barrel
  //2 forward cd
  //3 backward cd
  //-1 error

  //experimenttype 
  //0 one cd forward
  //1 one cd backward
  //2 no cd
  //3 barrel + 1 cd backward
  //4 barrel + 2 cds

  //barreltype
  //0 old no BE detectors
  //1 new

  //cdpad
  //0 no cd pads
  //1 cd pad detectors

  if(ExperimentType == 0){
    if(HCID < 4)
      return 2;    
  }
  if(ExperimentType == 1){
    if(HCID < 4)
      return 3;    
  }
  if(ExperimentType == 2){
    if( (HCID-FirstBarrelHC)%Quadrants == 0)
      return 0;
    else if ( (HCID-FirstBarrelHC)%Quadrants == 2)
      return 1;
    else{
      cerr<<"Error in HitCollection ID: " << HCID << endl;
      exit(5);
    }
  }
  if(ExperimentType == 3){
    if(HCID < FirstBarrelHC){ // cd 0-3 BCD 4-7 BCDE
      if(HCID < 4)
	return 3;
      else
	return -1;
    }
    else{
      if( (HCID-FirstBarrelHC)%Quadrants == 0)
	return 0;
      else if ( (HCID-FirstBarrelHC)%Quadrants == 2)
	return 1;
      else{
	cerr<<"Error in HitCollection ID: " << HCID << endl;
	exit(5);
      }
    }
  }
  if(ExperimentType == 4){
    if(HCID < FirstBarrelHC){ // cd 0-3 FCD 4-7 FCDE 8-11 BCD 12-15 BCDE
      if(CDPad){
	if( HCID%2 == 0) 
	  return HCID/8+2; // =2 FCD =3
	else
	  return -1;
      }
      else // cd 0-3 FCD 4-7 BCD
	return HCID/4+2;
    }
    else{
      if( (HCID-FirstBarrelHC)%Quadrants == 0)
	return 0;
      else if ( (HCID-FirstBarrelHC)%Quadrants == 2)
	return 1;
      else{
	cerr<<"Error in HitCollection ID: " << HCID << endl;
	exit(5);
      }
    }
  }
  return -1;
}
